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We investigate the asymptotic properties of inertial modes confined in a spherical shell 
when viscosity tends to zero. We first consider the mapping made by the characteristics 
of the hyperbolic equation (Poincare's equation) satisfied by inviscid solutions. Charac- 
teristics are straight lines in a meridional section of the shell, and the mapping shows 
that, generically, these lines converge towards a periodic orbit which acts like an attrac- 
tor (the associated Lyapunov exponent is always negative or zero). We show that these 
attractors exist in bands of frequencies the size of which decreases with the number of 
reflection points of the attractor. At the bounding frequencies the associated Lyapunov 
exponent is generically either zero or minus infinity. We further show that for a given 
frequency the number of coexisting attractors is finite. 

We then examine the relation between this characteristic path and eigensolutions of 
the inviscid problem and show that in a purely two-dimensional problem, convergence 
towards an attractor means that the associated velocity field is not square-integrable. We 
give arguments which generalize this result to three dimensions. Then, using a sphere 
immersed in a fluid filling the whole space, we study the critical latitude singularity and 
show that the velocity field diverges as l/Vd, d being the distance to the characteristic 
grazing the inner sphere. 

We then consider the viscous problem and show how viscosity transforms singularities 
into internal shear layers which in general betray an attractor expected at the eigenfre- 
quency of the mode. Investigating the structure of these shear layers, we find that they 
are nested layers, the thinnest and most internal layer scaling with i? 1 / 3 -scale, E being 
the Ekman number; for this latter layer, we give its analytical form and show its simi- 
larity to vertical |-shear layers of steady flows. Using an inertial wave packet traveling 
around an attractor, we give a lower bound on the thickness of shear layers and show 
how eigenfrequencies can be computed in principle. Finally, we show that as viscosity 
decreases, eigenfrequencies tend towards a set of values which is not dense in [0,2O], 
contrary to the case of the full sphere (f2 is the angular velocity of the system) . 

Hence, our geometrical approach opens the possibility of describing the eigenmodes 
and eigenvalues for astrophysical/geophysical Ekman numbers (10~ 10 — 1CP 20 ), which 
are out of reach numerically, and this for a wide class of containers. 



2 



M. Rieutord, B. Georgeot and L. Valdettaro 



1. Introduction 

Incrtial waves, which propagate in rotating fluids thanks to the restoring action of the 
Coriolis force, can generate very singular fluid flows when they are confined in a closed 
container. These very special properties of incrtial modes were first noticed in the theo- 
retical work of K. Ste wartson and others (fgtewartson fc Rickard 1969 ; [gtewartson 1971 



1972a| ,|E|; jWalton 1975|; London fc Shen 1979|). T hey appeared again recently in numerical 



investigations by HoUerbach & Kerswell (1995), Rieutord (1995), Rieutord & Valdettaro 
(1997f ), |Fothcringham fc HoUerbach (1998 ) and show an even greater generality since 



they are also present in stratified fluids (Maas & Lam 1995; Rieutord fc Noui 1999] ) or 
rotating stratified fluids (Dintrans et al. 1999). 

The particularity of all these waves (inertial, gravity, gravito-inertial) is that their 
associated modes are solutions of an ill-posed boundary-value problem when they are 
confined in a close container: the partial differential equation is of hyperbolic or mixed 
type in the spatial variables. This yields all kinds of singularities. When viscosity is 
included, these singularities are regularized but they still play a central role in featuring 
the shape of inertial modes of a rotating spherical shell: in particular, they control the 
asymptotic limit of small diffusivitics which is the relevant limit for astrophysical or 
geophysical applications. 

The aim of this paper is to present what we believe to be the asymptotic limit of 
inertial modes in a spherical shell when viscosity tends to zero. In the first part of the 
paper we shall present the main features of the solutions of this problem when viscosity 
is omitted. For this purpose we examine the trajectories of characteristics in a meridional 
plane of the shell as if they were trajectories of a dynamical system in some configuration 
space. We then focus on the relation between these trajectories and the eigenfunctions 
in two and three dimensions. We end this part with a close look at the critical latitude 
singularity. In the second part we investigate the changes brought on by viscosity and 
we examine more closely the structure of shear layers which arise. Then, by studying the 
behaviour of a wave-packet, we show how eigenvalues and eigenmodes may be computed 
in the asymptotic limit of a small viscosity. We conclude this part by a brief discussion of 
the distribution of eigenvalues in the complex plane. The paper ends with a discussion of 
the more general cases including containers with a different shape and of the applications 
of the present theoretical results. 

As this paper is rather long and goes through some mathematical developments which 
may be skipped at first reading, we suggest the casual reader to skip subsections §2.2.1- 
5, 2.3.1-2 and 2.4.1-2 and be lead by the introductions of sections 2.2, 2.3 and 2.4 and 
then jump to 2.5 and 2.6. The second part is not so mathematical but the details of the 
boundary layer analysis (§3.2.2-4) can be skipped at first reading. 



2. Some properties of inviscid solutions 

2.1. Equations of motion 

We consider a fluid with no viscosity contained in a spherical shell whose outer radius is R 
and inner radius r/i? with r/ < 1. The fluid is rotating around the z-axis with the angular 
velocity fL Using (2f2) -1 as the time scale and R as the length scale, small amplitude 
perturbations obey the linear equation 

§+e z xu 

V • u = 
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where u is the velocity held of the perturbations and p is the reduced pressure perturba- 
tion. The boundary conditions are simply 



u ■ e r 







at r — r\ and 



(2.2) 



As in Rieutord & Valdettaro (1997) and Dintrans, Rieutord & Valdettaro (1999) we shall 



e q will denote the 



use spherical coordinates (r, 6, <p) or cylindrical coordinates (s, ip, z 
unit vector associated with the coordinate q G {r, 6, ip, s, z}. 

Wh en t he time dependence of the solutions is assumed proportional to exp(icji), equa- 
tions (2.1) may be cast into a single equation for the pressure, namely 



vV 



1 d 2 p 

lu 2 dz 2 







(2.3) 



which has been referred to as Poincare equation since the work of Cartan (1922 ). This 
equation is completed by the boundary condition u ■ e r — which reads 



— uj 2 e r ■ Vp + iuj(e z x e r ) ■ Vp + (e z ■ e r )(e z ■ Vp) = (2.4) 

when expressed with the pressure; it applies at r = r\ and r = 1. 

As is well-known (Greenspan 1969), the Poincare equation is hyperbolic since for all 
modes lu < 1. Therefore, a first step in the analysis of this equation is the determination 
of the characteristics surfaces; for this purpose, we note that second-order derivatives of 
this equation read 

d 2 p d 2 p a 2 d 2 p 
dx 2 dy 2 co 2 dz 2 

where a 2 = 1 — to 2 , x and y are the cartesian coordinates in a plane z — 0. These 
second-order derivatives define characteristics surfaces such as z — f(x,y) = where / 
verifies : 



d£\ 2 (d£\ 2 = rf 



dx J ' \dy ) ui 2 ^ 
These surfaces are known as 'surfaces of constant slope' as any tangent plane makes the 
same angle $ = ^ — arcsinw with the equatorial plane z — 0. In fact, these surfaces may 
be generated by a family of such planes or by cones which aperture angle is 



A = arcsinoj (2-6) 

which is also the critical latitude in a sphere. In a meridional plane, the trace of these 
surfaces are simply straight lines making the angle A with the rotation axis. 

A second step in the analysis of the Poincare equation is to examine the separability of 
the variables. Because of the symmetry of the problem with respect to rotations around 
the 2-axis, the ip variable may always be separated from the two others. This implies 
that solutions may always be expressed as 

]>>„M)e"^ 

m 

and that each Fourier component p m (r, 9) is independent of the others. 

The two other coordinates, however, are not separable in the general case. This point 
may be understood easily if we recall that through a linear transformation of the z- 
coordinate, the Poincare equation may be transformed into the Laplace equation as first 
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Figure 1. (a) The four directions of propagation and the corresponding four values of S n 's. (b) 
A sketch of the notations used to describe the mapping. The numbers (1,4), (2,3) etc. indicate 
the possible directions of propagation of characteristics as shown in the panel (a). 



shown by Bryan 



(see also Greenspan 1969). In this transformation the boundaries 



transform into one-sheet hyperboloids. One therefore adopts an ellipsoidal coordinate 
system within which one of the boundaries is a surface of coordinate; unfortunately, 
since the two bounding hyperboloids are not confocal, coordinates can be only separated 
on one of the boundaries of the spherical shell. Hence, one may simplify, and actually 
solve, the problem either in the full sphere (see Greenspan 196S ) or in the infinite fluid 
outside a sphere (see below). 

We have now seen the basic ingredients which make this problem difficult: hyperbolicity 
(ill-posedness) and non-separability. 

2.2. Orbits of characteristics 

As shown by the foregoing discussion the difficulty of the problem lies in the behaviour 
of the solutions with respect to the coordinates in a meridional plane (s, z) or (r,0). We 
shall therefore restrict our analysis to this plane where characteristic surfaces are simply 
straight lines; thanks to the exp(im</?)-dependence, our results will apply equally to ax- 
isymmetric or non-axisymmetric modes. Indeed, using the separation of the ^-variable, 
we eliminate second-order derivatives in tp and characteristic surfaces are just cones in- 
dependent of m. We shall therefore study, in the fo llowing subsections, th e trajectories 
of characteristics in a meridional plane as we did in Dintrans et al. (1999 ). Since this is 
a rather technical matter, the reader may first skip it and directly jump to section 2.3 
where the results are summarized. 

2.2.1. The mapping 



From (2.5) we derive the well-known equations of the two families of characteristics: 



ujz ± as = (2-7) 

where u± will designate the characteristics coordinates. u+ and u- are constant along 
characteristics of positive and negative slopes, respectively. 

To describe the paths followed by characteristics, we first study the map which relates 
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Figure 2. The resulting mapping in the case of a shell with r\ — 0.35 when the frequency 
is lu — 0.40782. Twelve points of discontinuity indicate the projection of the 'shadow' and the 
critical latitude of the inner sphere on the outer one (see text); the four apparent discontinuities 
due to the periodicity in the [0,47r] interval are not counted. 



the position of the n+l th reflection point to the one of the n th both taken on the outer 
sphere. For this purpose we mark out these points by an angle (f> G [0, 2ir] which is 
identical to the latitude when < <f> < 7r/2; this is indeed more convenient than the 
colatitude. If one reflection is needed on the inner sphere then we have 



(2.8) 



sin(</>„ - 5 n ) = r?sin(/3 - S n ) 
sin(<^ n+1 - 5 n+1 ) =r/sm(P - 5 n+1 ) 

where j3 is the position of the reflection point on the inner shell and 8 n = ±7r/2 ± A is 
the direction of the characteristic which can take four values as illustrated in figure [j]a. 
If the two reflection points are simply connected by one segment of characteristics then 
the recurrence relation is either 

<t> m + <t> m +i = -2A [2tt] (2.9) 
when the characteristic has a positive slope, or 

= 2A [2tt] (2.10) 
when the characteristic has a negative slope. These notations are summarized in figure 0b. 



From the expression (2.8), (2.£) and (2.10) one can compute the map 



'M+i — 



/+(0n) 

f-(<M 



Such a map is bi-valued since one may compute the image by first applying a positive- 
(/+) or a negative- (/_) slope characteristic. However, such a representation is not con- 
venient for iterating the map since for each iteration one has to decide whether to use 
/+ or /_ . We therefore constructed a single- valued map which is defined in the following 
way: Considering a point of the outer sphere, we mark it out by the angle <p € [0, 27r] if 
it is to be iterated by /_ or by the same angle plus 2tt if it is to be iterated by /+. We 
thus define the map: 
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Figure 3. The third (left panel) and eighth (right panel) iterate of the map in the case of a 
mode with lo — 0.40782 for a shell of aspect ratio r\ = 0.35: note the appearance of fixed points 
marked out by the intersection of the curve with the line <f>„+ P = 4>n ■ They indicate the existence 
of attractive periodic orbits of period 3 and 8 respectively. Note also the increased number of 
discontinuities of the mapping. 



/ : [0,4tt] — ► [0,4tt] 

4>n * f{4>n) = 4>n+l ( 2 -H) 

which is one-to-one except at some points of discontinuity and which can be easily iter- 
ated. An example of this map is given in figure |^. 

One of the remarkable features of this map is that it is not continuous. The disconti- 
nuities occur at the colatitudes (in the first quadrant) 

6± = A ± arcsin.77, 6 S = A + arcsin(?7 cos 2 A) 

9± are delineating the 'shadow' projection of the inner shell on the outer shell (see 
figure [5] for an illustration of the shadow) . They illustrate the case when a characteristic 
is tangent to the inner sphere at critical latitude. 9 S is the colatitude of the projection 
of the inner sphere's critical latitude on the outer sphere. These angles (6±) delimit the 
intervals where the map is contracting |/'| < 1 or dilating |/'| > 1 or neutral |/'| = 1. 
When the map (2.11) is iterated as in figure || and since, generically, discontinuities 



are not mapped onto themselves[f], their number increases proportionally to the number 
of iterations; also some fixed points appear which indicate the existence of attractors, i.e. 
attractive periodic orbits which we discuss below (§ [2. 2. 4 ). We would therefore expect that 



the basins of attraction of the infinitely iterated map, containing an infinite number of 
intervals at smaller and smaller scales, would have a fractal structure; however, numerical 
studies show only isolated accumulation points which are actually the fixed repulsive 
points of the mapping, i.e. the repellors (see figure ||). 

We therefore see that the structure of basins of attraction is much more complicated 
than in the case studied by Maas fc Lam (1995| ) and may represent the general case for 
such systems. 



f The case when all dis continuities are mapped onto themselves corresponds to periodic orbits 
of the shadow (see § |2.2.3 ). 
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Figure 4. Zoom of the N th iterate (N = 1200) of the map_near an accumulation point. The 
aspect ratio n and the frequency uj are the same as in figure |3[ On the left panel we clearly see 
the presence of the accumulation point 4> R ; the two arrows indicate the basins of attraction of 
the polar attractor in figure 0a; the other segments belong to the basin of the other (equatorial) 
attractor. On the right panel we see that the width of the intervals of the basins of attraction 
vanishes geometrically as the accumulation point is approached. 



2.2.2. Orbits and Lyapunov exponents: the full sphere case 

Once the mapping is known, we may examine the trajectories of characteristics. For 
the sake of clarity, it is useful to first consider the case of the full sphere for which only 
( [2.9[ ) and ( 2.10 ) are necessary. 

We first note that the number of reflection points of a periodic orbit is necessarily 
even, i.e. 2q, when the starting point is not a critical latitude. Applying alternatively 
(pjj) and Q2.10D, we have 



h = -02 + 2A [2tt] 
-(f> 2 = 03 + 2A [2rr] 

hq-l = -4>2q + 2A [2tt] 

-02, = 0i + 2A [2tt] 



(2.12) 



Summing all these equations leads to the conclusion that a periodic orbit is such that 



_ pn 



(2.13) 



where p and 2q are integers which represent the number of crossings of the orbit with 
respectively the axis of rotation or the equator. 

From the preceding results, it turns out that all orbits such that A = m with r irra- 
tional, are ergodic (quasi-periodic) . At this point it is worth noting that eigenfrequencies 
of inertial modes in the full sphere are in general associated with quasi-periodic orbits. We 
prove in appendix ^ that, for instance, the first axisymmetric inertial mode of frequency 
\/3/7 is associated with a quasi-periodic orbit. 

To co nclud e with the full sphere, we just need to point out that thanks to relations 
( [2.9| ) or ( 2.1 0| ) , it is clear that the Lyapunov exponent is always zero. Indeed, if we recall 
the definition of the Lyapunov exponent A associated with an orbit: 



A 
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Figure 5. The shadow pattern for a spherical shell with 77 = 0.2 when A = 7r/6 (u = 0.5). 



2.2.3. Orbits and Lyapunov exponents: the shell case 

We now turn to the shell case. As we already observed the map has discontinuities 
defined by the shadow of the inner sphere on the outer shell and the projection of critical 
latitudes. If the inner sphere is sufficiently small, the shadow may draw a periodic pattern 
if the critical latitude A is commensurable with it. A simple case is illustrated in figure ||. 
For these frequencies, orbits starting in the shadow remain in the shadow while those 
starting outside remain outside. In this way, one can construct the set of all periodic 
orbits with A = i.e. all neutral periodic orbits. 

The fact that periodic orbits outside the shadow are neutral is obvious; the case of 
orbits in the shadow is less obvious but we may consider the fact that if these orbits 
were not neutral, the shadow would not map onto itself. As a short exercise, we may 
follow one such orbit for u> — sin (71-/ 2 (2p +1)). The shadow should cross only twice the 
location of the inner sphere. It bounces first on the inner sphere, at an angle ft, and 
then on the outer sphere at the angle ft given by sin(ft — 5) = 77 sin(ft — 8) (5 is one 
of the angles S n in figure Then it bounces 2p times on the outer sphere. After each 
two reflections, the angle (j> changes into 4> + AS. Therefore after 2p bounces, the angle 
on the outer sphere is ft + ApS. It bounces then again on the inner sphere, hitting it at 
the angle ft such that sin(ft + S + 4p5) = 77sin(ft + S). But 2(2p + 1)5 = n therefore 
6+4pS = n-6[2ir]. Sosin(ft-<5) = -yysiftft+J) = T)sin(0i-6). Therefore ft = -ft[27r] 
or ft = 7T + ft + 2<5[27r] . Repeating this entire process 2(2p + I) times, one comes back 
to the original ft . 

In fact periodic orbits of the shadow do not exist for all A's commensurable with tt. 
Indeed the image of the shadow must not be split by discontinuities; this implies that 
p or q cannot be too large or the shell too thin. More precisely, for a given 77, periodic 
orbits of the shadow exist if: 



it is clear that for the full sphere 



dct> n +i 



= 1 Vn, so A = 0. 



arcsin?/ < A < arccosry (2-14) 
If 77 > 1/V2, only one neutral periodic orbit exists: it is such that A = 7r/4or uj = 1/V2. 



Inertial waves in a rotating spherical shell: attractors and asymptotic spectrum 9 





iprrj 



3. -0.5 - 




arcsin(u)/TT 



Figure 6. The lyapunov exponent of one basin as a function of frequency (left) or critical 
latitude (right). Note the symmetry of this last plot with respect to n/4. The vertical dotted 
lines mark the critical latitudes it/6, 7r/4 and 7r/3; the aspect ratio of the shell is n — 0.35. 



More generally, for a given aspect ratio rj, one may determine all the frequencies associ- 
ated with neutral periodic orbits and their number increases as the size of the inner shell 
decreases. These frequencies are obviously determined by ( 2.13 ) but due to the finite size 
of the shadow, one needs to eliminate large values of p and q. If we remark that the most 
robust periodic or bits ( when rj increases) are those with p = 1 (restricting ourselves to 
A < 7r/4), relation ( 2.14| ) easily bounds the permitted values of q. For a shell like the core 
of the Earth, where rj = 0.35, only q = 2, 3, 4 are possible. 

The frequencies of neutral periodic orbits are important as they shape the form of the 
Lyapunov exponent curve since, then, intervals contracting through /+ are exactly dilated 
by /_ which makes the Lyapunov exponent vanish. As a consequence, frequencies in the 
neighbourhood of one such frequency are associated with very long orbits having a small 
Lyapunov exponent. This is the reason why 'windows' appear near these frequencies, 
especially near lu — l/y/2 (A = 7r/4) as clearly shown in figure [| 

This latter figure shows many spikes which in fact betray the presence of other periodic 
orbits which we shall call attractors after Maas & Lam (1995). For such orbits A < 0. 
In fact for this system, all the orbits (except may be some isolated ones) verify this 
inequality and no chaos is possible: the configuration space being one-dimensional and 
the mapping being one-to-one. 

2.2.4. Some properties of attractors in the shell 

In order to make the dynamics of attractors clearer, it is useful to concentrate on a 
specific example which can be computed explicitly. For this purpose, we choose a spherical 
shell similar to that of the core of the Earth for which r] = 0.35 and we focus on the 
orbit with four reflections on the outer or inner shell as shown in figure |. If the shell is 
thinner, this orbit is localized in the vicinity of the equator which is the kind studied by 



Btcwartson (1971| , |1972a 

The Lyapunov exponent of this orbit may be computed explicitly in the following way: 
Let us first recall that plane inertial waves reflecting on a surface oriented by ft verifies 
the relation 
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Figure 7. (a) Two attractors coexisting when ui = 0.4051; we call them 'polar' and 'equatorial' 
respectively as they are characterized by the fact that the reflection on the inner sphere occurs 
above or below the critical latitude; the arrows indicate the direction of focusing, (b) Lyapunov 
exponent as a function of frequency in the vicinity of to — 0.4051; the thick solid line denotes the 
equatorial attractor while the thick dashed line rep rese nts th e po lar attractor; these two thick 
lines have been derived from the analytic formulae ( B 4 ) and ( B 8 ) given in appendix [b| 



where fcj and k r are the wave vectors of the incident and reflected waves respectively 
( preenspan 196S ). When applied to a reflection on a sphere, this relation implies 

sim>± A) 

K T — Ki . . . 

sm(0 =p A) 

where <fi is the 'latitude' of the reflection point. From this relation, we can compute the 
rate of contraction of an infinitesimal interval in the neighbourhood of a reflection point 
of the orbit. Therefore, the Lyapunov exponent of an attractor with N reflection points 
is simply given by: 



N 

A 



N sm(& t A) v I 

where the cafe's describe the reflection points of the attractor (a periodic orbit). At this 
point we shall define two useful quantities characterizing attractors, namely their length 
and their period. We shall call the number of reflection points (on both spheres) the 
length of the attractor while its period will refer to the minimum number of iteration of 
the mapping needed to generate its associated fixed points; because of the definition of 
the mapping, the period is also the number of reflection points on the outer sphere. For 
instance, the period of the equatorial attractor of figure [?]a is '3' and its length is '4', 
while for the polar one, the length is TO' and the period '8' (we do not allow reflections 
on the axis). We give in appendix |b] the explicit calculations relevant to the orbits of 
figure ^ 

The curves A(a;) (figure 0b) show the same feature: in the interval of existence of the 
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orbit, A varies between and — oo. The two extremes correspond respectively to the cases 
when the reflection on the inner shell occurs at the equator or at its critical latitude. We 
note that the vanishing value of the Lyapunov exponent does not mean that for this 
frequency the orbit is no longer an attractor: it simply means that convergence towards 
the attractor is no longer exponential; in most cases, it does converge, but algebraically. 

For later use, we also computed the behaviour of A(a>) in the vicinity of the point loq 
such that A(ujq) — 0. We find in the particular case of the equatorial attractor that 



A(lu) = -K(lu -u ) 1/2 , with K = 4.8184, lu = 0.403112887 (2.16) 

In fact this behaviour is general as is shown in appendix |b[ Let us also emphasize that 
when A = 0, the orbit is just one broken line connecting: 

• the equator or the pole of the spheres to another equator or pole, 

• the equator or the pole of the spheres to the critical latitudes of the outer sphere. 
With the terminology in use for dynamical systems, such an orbit restricted to the first 
quadrant is self-retracing: a mass-point would go back and forth on the same trajectory. 

The foregoing results make the shape of figure ^ quite clear now. Most of the 'spikes' 
shown in this graph will therefore tend to — oo as the number of points of the graph is 
increased. But to be complete, we need also mention some cases when a segment of an 
orbit intercepts the inner shell after being tangential to it. In this case the curve A(w) 
has a discontinuity and does not reach — oo. 

We therefore see that attractors are featuring figure ^. As it will be clear later they 
also feature the shape of the asymptotic spectrum. It is thus interesting to know some 
elementary properties of these geometrical objects. 



From a rather large number of computations we observed, as [Maas fc Lam (1995|) , that, 
in the first quadrant, not more than two attractors may coexist for a single frequency. 
However, these two attractors can be used to construct other attractors which are just 
their image symmetrized with respect to axis of rotation and equator. Considering the 
propagation of characteristics in the full meridional section of the shell, we observe that 
these lines can converge towards six attractors at most. Using the properties of the 
mapping ( [2.11 ), we have been able to prove under certain hypotheses (see appendix ^|), 



that the number of attractors is bounded by the number of points of discontinuity (which 
is twelve). 

We also computed the interval of existence, in frequency space, of a large number 
of attractors so as to show its relation with the length of the attractor. As shown in 
figure |8j the interval of existence is well correlated with the inverse square of the length. 
We explain this law in the following way: for a very long attractor of length N 1 
the number of reflections on the inner and outer shell scales with N, therefore the mean 
angular distance between the critical latitude and the nearest reflection point is 0(1/N). 
This implies that just a 0(1/N 2 ) variation in frequency is necessary to shift this point 
to the critical latitude. 

The latter result has an interesti ng co nsequence on the shape of the curve A(w) for a 
given long attractor. Indeed, from ( PU51 ) the diver gence toward — oo of an attractor of 
length N is of the form: 

A~iln[tf(w-w c )] 

where lu c is the frequency of the singularity of A; we used the fact that ( |2.15| ) is dominated 
by one term ( <pk — A ~ with cf>k — A ~ N(u> — u c ) ); since the interval of existence of the 
attractor scales like 1/N 2 , if we choose a point such that uj — uj c = a/N 2 the Lyapunov 
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Figure 8. Interval of existence in latitude 
of several different attractors, plotted as 
a function of their length N The straight 
solid line is the inverse square of the length. 



Figure 9. Distance to the point on the ex- 
ternal sphere at critical latitude, for several 
different attractors, plotted as a function of 
the length N. The straight solid line is the 
inverse square of the length. The dashed 
line gives the lower bound (in distance) for 
physical ly r elevant attractors (see discus- 
sion in ^3.1). 



exponent will scale like A ~ — ln (^v) w hlch vanishes at large N. This means that the 
singularity of the Lypaunov curve occupies a smaller and smaller fraction of the interval 
of existence of the attractor. Hence for long attractors, the Lyapunov exponent is very 
small in a larger and larger part of their interval of existence. This explains why long 
attractors appear numerically as weakly attracting eventhough their Lyapunov exponent 
may diverge. 

In figure ^ we show the distance of attractors to the point on the external sphere at 
critical latitude as a function of the length. Since a periodic orbit exists in a range of 
frequencies (see figure |J), instead of showing a single point we represent a vertical segment 
connecting the minimum and maximum distance over the entire range of existence of the 
attractor in frequency space. We see that the maximum distance is well correlated with 
the inverse square of the length. This distance is important in the final appearance of the 
attractor when viscosity is included. As it will be shown later on an example, the build- 
up of energy along an attractor can be impeded by the boundaries; this effect therefore 
puts an upper bound on the viscosity for the attractor to be visible. The dashed line 
in the figure gives the lower bound (in distance) for physically relevant attractors (see 
discussion and the example in §3.1). 



2.2.5. Differences with billiards 

It is interesting to compare the mapping defined by (2.8-2.10) with the billiard problem 
studied in classical chaos, where a particle bounces specularly on the walls of a cavity. 
Billiard phase spaces are two-dimensional (position and velocity direction) while the 
phase space of our problem is one-dimensional (in projection), since the only variable 
is the position along the circles representing inner and outer shells. The problem is not 
Hamiltonian, there is no conservation of the symplectic measure in phase space, and 
attractors and repellors exist. 

For the full sphere, rj — 0, we have seen that all the orbits are neutral (A = 0) and 
are either quasiperiodic (and ergodic) or periodic (when A is a rational fraction of tt). 
When rj is made nonzero, all quasiperiodic orbits are instantaneously destroyed, and the 
periodic orbits remain neutral until they are eventually destroyed when rj increases, the 
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rational values with smaller denominator surviving last. It is interesting to note that 
this situation is exactly the opposite of the Kolmogorov-Arnold-Moscr (KAM) theorem 
valid for Hamiltonian systems close to integrability (see Arnold 1989). In this latter 
case, it is well-known that for an integrable Hamiltonian system all orbits lie on tori, and 
orbits are organized in families, either quasiperiodic (and ergodic) or periodic, for a given 
torus. If one perturbs such an integrable Hamiltonian system by a sufficiently smooth 
perturbation, the KAM theorem states first that all rational tori (with periodic orbits) 
are instantaneously destroyed, and second that irrational tori (with quasiperiodic orbits) 
disappear one after the other when the perturbation is increased, the last to disappear 
being the 'furthest' to the rationals. 

2.3. Relations between orbits of characteristics and eigenf unctions 

In the preceding subsections we have shown that, generically, characteristics converge to- 
wards attractors which are formed by a periodic orbit. These attractors live in a frequency 
band whose size decreases with the length (or period) of the attractor. The attracting 
power is measured by a negative Lyapunov exponent which generically varies between 
and — oo when the frequency band is scanned. Several (less than 6) attractors may 
coexist for a given frequency; in this case they own each a basin of attraction (described 
for instance as a set of points on the outer boundary) whose structure is governed by 
accumulation points (see fig. ||). 

We have also found periodic orbits which are not attractors; their frequency can be 
written u) = sin(pir/2q) where (p, q) are chosen in a finite set of integers. The number of 
these orbits is therefore finite but increases as the radius of the inner core decreases. The 
frequencies of these orbits will prove to be useful since in their neighbourhood, attractors 
have very great length and, therefore, very small (in absolute value) Lyapunov exponents. 
This will influence the shape of the asymptotic spectrum (i.e. with low viscosities). 

Now we shall see how the eigenfunctions are influenced by the presence of an attractor. 

2.3.1. The two-dimensional case 

Because of the simple form of the Poincare equation in two dimensions, which may be 
written as 



d 2 p 

du+du- 



(2.17) 



early investigations have focused on this case in particular those of mathematicians 
The relevant contributions are those o f [Bourgin fe Puffin (1939| ), |john (194l| ), |H0iland| 



(1962|) , [Franklin (1972| ) , |Ralston (1973[ ) and ^chacffcr (1975|). M uch of this previous work 



is concentrated in a theorem demonstrated in Bchacffcr (1975 ) which states that: 

There are non-trivial solutions of ( 2.1 r \ j if and only if there exists an integer n such 
that all reflected rays close after precisely 2n reflections. If there is one solution, then 
there are infinitely many, linearly independent solutions. 

In other words, eigenvalues of regular modes are always associated with periodic orbits 
and these eigenvalues are always infinitely degenerate. Since the reflected rays must close, 
starting from any point of the curve, the Lyapunov exponent of such periodic orbits is 
always zero. Note the difference with the three-dimensional case where eigenmodes in the 
full sphere are associated with ergodic (quasi-periodic) orbits (cf § |2.2.2 and appendix |X]) 
and eigenvalues are non-degenerate. 

Another interesting result was derived from mathematical analysis by Ralston (1973| ) . 
Namely, it states that the velocity field associated with a solution of ( |2.17 ) is not square- 
integrable when characteristics are focused towards a wedge formed by the boundaries. 
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An example of such a singular flow is given in Wunsch (1968) with the case of internal 
waves focused by a sloping boundary. In the interval of frequencies where the velocity 
field is not square-integrable, eigenvalues do not exist and the point spectrum of Poincare 
operator is said to be empty. 



The foregoing results may be generalized to our case, or the one studied by Maas & Lam 



(1995), where characteristics are focused towards an attractor. Indeed, let us consider 
the total kinetic energy of a 'mode' associated with an attractor; using characteristic 
coordinates, this quantity reads 

1=1 \\v\\ 2 du + du- + / ||w|j 2 iiu+du_ 
Js Jcs 

where S designates a neighbourhood of the attractor and CS the remaining 'volume'; we 
assume that the limits of S are made up of characteristics. If the attractor is of length 
N then this integral can be split into N pieces 

JV 



/ = J \\v\\ 2 du + du_ 



where we neglected the contribution from CS. Now each of these pieces can be split again 
into an infinite number of rectangles Rk with sides made up of characteristics. Hence we 
write 



||w|| 2 d?i+(iu_ = / 1 1 v\ | 2 du + du _= Ik 

k=l Rk k=l 

In the vicinity of the attractor, these rectangles are very elongated: one side remains 
0(1) long while the other shrinks to zero as the attractor is approached. 

Now we take the two long sides as images through the mapping made by the charac- 
teristics. The mapping has a contracting rate given by e A < 1 where A is its Lyapunov 



exponent. Maas & Lam (1995) have shown how one can construct the stream function 
in the whole domain by iterating an arbitrary function given on its boundary. When 
the attractor is approached, the scale of the stream function vanishes while its ampli- 
tude remains constant; therefore the kinetic energy is amplified by a factor e~ 2A at each 
iteration of the mapping. Noting that one rectangle is smaller by a factor e A than its 
predecessor, we can derive the iteration rule 

4+1 = e- A I k 

which shows that the integral I is infinite. We may note in passing that if dk is the 
distance of the k th characteristic to the limit cycle, then dk = doe kA while the amplitude 
of the velocity field is Vk = voe~ kA . This shows that the velocity field diverges as the 
inverse of the distance to the attractor. 

Therefore, as in the case of a wedge, the velocity field is not square-integrable when 
characteristics converge towards an attractor. 

2.3.2. The three-dimensional case 

The 3-D case has not benefitted from the same interest by mathematicians. In this case 
the Poincare equation contains first or zeroth order derivatives which cannot be elimi- 
nated. Let us rewrite it using cylindrical coordinates and assume a exp(imip) dependence 
of the pressure; thus 

d 2 P 1 dP a 2 d 2 P m 2 P 



ds 2 s ds ui 2 dz 2 



= 
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Figure 10. A sketch for the illustration of Riemann's method; the lines (SP) and (SQ) are 

segments of characteristics. 



The canonical form of this equation is obtained using characteristics coordinates: 



o 2 p 

du+du- 



dP dP 



du+ du- 



-P = 



(2.18) 



which is known as the Euler-Darboux equation ( Dautray fc Lions 1984-1985 ) 



A general solution of Euler-Darboux equation may be obtained with Riemann's method 
( polombo 1976 ; Zwillinger 1992). With this method one may express the value of the 
solution at one point when 'initial' data are given on an arc joining two points on charac- 
teristics 'emitted' from the point considered (see hgure 10). However, one needs to know 
the Riemann function (which plays an equivalent role to the Green's function of elliptic 
problems). To determine this function it is useful to rewrite the pressure fluctuation as 
P = IL/y/s; doing so, the first derivatives of Poincare equation are eliminated but are 
replaced by the term !I/4s 2 . The equation for II is therefore: 



d 2 U 



n 



= o 



The Riemann function is a solution of this equation^ which meets the additional con- 
ditions: 



R{u+,u) = 1, 



dR 

du + 



or, equivalcntly, 



dR 

du- 



= 



(2.20) 



R(u' + ,u-) = R{u +1 u'_) = 1 V{u' + ,u'_) E V 



f In fact, it is a solution of the corresponding adjoint operator which is identical to the original 
in this case. 
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We shall call S(u+, u_) the point where the solution is computed and M(u', , u'_) a point 
running on the arc of data; T> is the area defined by (SPQ). 



As Friedlander & Hcins (1968) have noted, (2.19) is invariant for all the transformations 



leaving 

(u' + — u + )(u'_ — U-) 
(u' + — u'_)(u + — U-) 



z 



invariant. Therefore, seeking a solution of the form R(z), one finds that this function 
verifies the differential equation: 

z(l - z)R" - (2z - 1)R' + fi(/j, -1)R=0 

where we set [i = m + 1/2. This is a special case of the differential equation of Gauss 
hypergeometric function, i.e. F([a, 1 — fi; 1; z); in fact, it is just the equation satisfied by 
Legendre functions of index — 1. 

Since z = z(S,M), we shall write Riemann's function as R(S;M). Hence the formal 
solution of the problem is 

n(S) = i(n(P)+n«3)) 

+ l Jpq R ^ M > - £:*-) + "< M > - SET*-) < 2 ' 21 ' 

which shows how the value of II at S may be constructed from the data given on the arc 
PQ. 

A simpler formula can be obtained for axisymmetric modes when the meridional stream 
function ip is considered^. After a similar transformation, where we set tp = y/s^!, we 
find that the associated Riemann function is also a Legendre function with /i = —1/2. 
If the arc PQ is taken on the boundary then W = and the expression ( |2.2l| ) simplifies 
into 

If / 9* \ 

*(5) = - / R(S; M) ^—du+ - ^—du- (2.22) 



2 J pq V du+ du- 



Let us now suppose that S is also on the boundary (on the inner sphere for instance) ; 
then = 0. If we introduce dC(*) = (j^du + - ™-duJ), then we have 



R(S;M)dC(ty) = 
PQ 



This relation holds also for neighbouring points (P'Q 1 S') of (PQS), thus 

f R(S';M)dC(V) = Q 

J P'Q' 

By subtracting these two equations we get 



R(S; M)dC(V) + / R(S;M)dC(V) + d<t>s 

PP' JQQ' JPQ 

f This function is such that u s = § §7 , u z — — 1 1^ 



f ^.(S;M)dC(*)=0 (2.23) 
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where dtps denotes the variation of the position of S. Since P' and Q' are in a neigh- 
bourhood of P and Q respectively and that R(S; P) = R(S; Q) = 1, the first integrals of 



(2.23) can be simplified so that: 

C(*)(P)+C(*)(Q) + ^f / §f(S;Af)dC(*) = (2.24) 

# JPQ V<P 

where Cf*l - d^_du+ _ d^_du^ 

Now we consider that PSQ are part of a limit cycle like the one of figure 0a (right); 
let us call T the fourth point of this cycle and let P n , S n ,Q n ,T n be the suite of points 



converging towards PSQT (i.e. those points at (03, $4, (j>i, ^2) respectively). ( |2.24j ) can 
be applied to the triangles P n , £„, Q ra and Q n , T n , P n +i and we get: 

C(*)(P n )=C(*)(P n+1 )+^ [ ™(T n ,M)dCm+^ I d -^{S ni M)dCm 

d( P JQnPn+1 °9 d( P JQnPn °<P 

The two integrals in the RHS are of order unity and we surmise that they do not 
cancel. Therefore the suite of C(^)(P n ) is diverging which means that the velocity fields 
tends to infinity when a limit cycle is approached. 

This result can be generalized to any limit cycle. We may also use it to generalize 
Ralston's theorem. However, in this latter case, it is more straightforward to note that if 
we are considering characteristics converging towards a wedge (for a three-dimensional 
problem), then in a neighbourhood of the apex of the wedge, first order derivatives are 



negligible compared to second order derivatives; therefore, in this neighbourhood, (2.18) 



can be transformed into (2.17) and Ralston's theorem applies. 



2.4. The critical latitude singularity 

The preceding discussion has shown (and proved in 2D) that the velocity field of 'modesfj 
associated with an attractor is not square-integrable: it diverges as the inverse of the 
distance to the attractor. 

We shall see now that this singularity is not the only one and that a milder one develops 
around the critical latitude of the inner sphere. |Stewartson fc Rickard (1969 ) were the 



first to notice that singularity and showed that it is integrable. Although the work of 
Stcwartson fc Rickard (1969| ) was restricted to the thin shell limit and was based on the 



use of Longuet-Higgins solutions of the Laplace tidal equation, we shall show that their 
result is in fact general. 



For quick reading demonstrations in §2.4.1 and §2.4.2 can be skipped 



2.4.1. The singular surfaces 

Let us consider a sphere immersed in a rotating fluid filling the whole space. We 
examine the oscillations of the fluid. Such modes are the corresponding modes of the full 
sphere when solutions regular at the origin are replaced by solutions regular at infinity. 

As for the full sphere, we use Bryan's transformation to convert the Poincare equation 
into the Laplace equation. Thus we set 

/ 

z = —t—z 
a 

To solve Laplace's equation, we therefore use ellipsoidal coordinates (£, 9, <p) similar to 
the spherical coordinates (for the angular variables 8 and <p): 

f We use quotes because modes refer usually to regular solutions with square-integrable ve- 
locity fields. 
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x = a cosh £ sin cos ip 

y = a cosh £ sin 6 sin <p )■ (2.25) 
z' = a sinh £ cos 

where a is the focal distance of the meridional ellipse. If we take the radius of the sphere 
as the unit of length, then a — I /a. We recall that using these coordinates, Laplace 
equation of an axisymmetric held transforms into : 

v2v = S5?| (" 8h{ f ) + SiW^ " (2 ' 26) 

whose solution, regular at infinity, reads: 



Qe(i sinh £)P e (cos 9) (2.27) 

where Qi is the second-kind Legendre function. 

If we now come back to the original coordinates, we may write the cylindrical coordi- 
nates (s, z) as 



s = — cosh £ sin ( 
a 



sinh £ cos 9 



and following the idea of |Greenspan (196E ) , we introduce 



(2.28) 



cost 



and i] = i sinh £ 



so that 



(2.29) 



as = ^(l-^)(l- v 2 ) 
luz = fir] 

The Jacobian of this new coordinate transform is 



(2.30) 



where we used 



J 



(2.31) 



c)fi 



A ' 



djt 

dz 



A 



dr] 

ds ~ 



a sr/ 



<>>! 

<)z 



A 



(2.32) 



with A = rj 2 — fj 2 . Hence, the transform is singular on the surfaces such that rj = ±/i. 
Since the solution Qi(r])Pi(fi) is regular in the fluid's domain, the singularity of the 
transformation makes the solutions singular when the coordinates map the space in a 
regular way. 

To discover which kind of surfaces hinds behind this equation (rj = ±/i), it is convenient 
to express £ (or rj) as a function of the cylindrical coordinates s and z. Eliminating 9 
from (|2.28|) and setting X = sinh 2 £, we find that 



X 2 + (1 + lu 2 z 2 - a 2 s 2 )X + lu 2 z 2 = 
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Figure 11. Meridional cross-section of the surfaces (cones) where the coordinate transform is 
singular. The signs +/— refer to the sign of the discriminant D. 



The solution of this equation gives the reciprocal transformation of coordinates (2.30) 
or ( 2.28D - When the roots are multiple, the transformation is singular; this happens when 
the discriminant D vanishes which is when 



D = (1 — ujz — as)(l — ujz + as)(l + ujz + as)(l + luz — as) = (2.33) 

One may easily verify that this equation is equivalent to rj 2 = fi 2 . 

We have therefore shown that the transformation is singular on four surfaces which 
are cones tangent to the sphere at the critical latitudes. This result is summarized in 
figure [ll]. 

2.4.2. The velocity field near the critical latitudes 

In order to present in a simple way the singularity of the velocity field near the 
critical latitude, we specialize our reasoning to the case of the tangent characteristic 
ujz = 1 — as which touches the sphere at s = a and z = ui. The velocity component 
parallel to this characteristic is such that Vji = uiv s — av z , or 

_ui 2 dP a dP 
a 2 ds ui dz 



Using (2.30) and (2.32), we find that 



a d[i 

y/l - T] 2 OP 



-(^Wl - + aWl - V 2 ) 

a or] 

Therefore, it turns out that if the right-hand side of this equation remains finite on 
the singular surface, then the velocity component Vji diverges as 1/A = \j\fT). In the 
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neighbourhood of the singular surface, D vanishes linearly with the distance to this 
surface; thus we see that the velocity field will diverge as one over square root of the 



distance to these surfaces as actually found by Stewartson & Rickard (1969) in the case 
of a thin shell. 

Let us show that the RHS of the latter equation is indeed finite. The characteristic 
ojz = 1 — as is such that fi — r\, therefore 

(l-/! 2 ) (dP dP" 



RHS = n 

a \o/j, on 

but — |^ = P^(rj)Qi(rj) — Pi(T))Q' e (rj) which is the wronskian of the Lcgendre 

functions; it is nonzero as Pi and Qg are linearly independent. 

Finally, using the same kind of arguments one may also prove that the component 
of the velocity field in the azimuthal direction is also singular while the component 
perpendicular to the singular surface remains finite. 

We therefore see that the velocity field possesses an integrable singularity but is not 
square-integrable; thus, if strictly quasi-periodic trajectories of characteristics exist, they 
would inevitably touch the critical latitude and their associated eigenfunction would be 
singular. Therefore no eigenvalue can be associated with ergodic trajectories of charac- 
teristics in a spherical shell. 

2.5. The toroidal (regular) solutions 

The foregoing two sections have shown us that inertial 'modes' of a sphericl shell hardly 
escape to singularities: one question therefore raises up: do regular modes exist at all? the 
answer is yes, indeed some regular solutions exist in the form of purely toroidal velocity 
fields associated with eigenfrequencie s w — l/(m + 1), m £ IN*. 



We pointed out these solutions in Rieutord fc Valdettaro (1997 ) but they appeared 



independently several times in the literature: Malkus (1967 ) noticed them while inves- 
tigating hydromagnetic planetary waves and Papaloizou fc Pringle (1978 ) called them 
'r-modes' because of their similarity with Rossby waves. 

However, the existence of these solutions is somehow puzzling since a plot of the 
trajectories of characteristics associated with these eigenvalues shows that most of them 
converge towards an attractor (for instance, when m = 2); how can we reconcile these 
two apparently contradictory facts? 

The answer lies in the specific form of the Poincare equation in these cases. Indeed, 
these modes are purely toroidal which means that for all points in the spherical shell 
we have e r ■ v = 0. From the expression of the velocity components as a function of the 



pressure fluctuation (see for instance Rieutord & Noui (1999)), this constraint can be 
transformed into the following equation 

uj dP mP z dP . . . 

—~s— 1 5 — = V(s, z) m the domain 

or as or ui Oz 

When this equation is combined with the Poincare equation, it turns out that the pressure 
must satisfy 

a 2 ( z d 2 P d 2 P\ m dP m 2 P 

■flfl-- ITT A ( 2 ' 34 ) 

osoz oz z ) lus os s 

The characteristics of this hyperbolic equation obey the differential equation: 



zdzds + s(ds) 2 = 
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They are therefore either straight lines parallel to the rotation axis ds = or circles 
parallel to the boundaries s 2 + z 2 = K. They cannot form orbits by reflections on the 
boundaries and therefore they do not impose any constraint on the solution; regular so- 
lutions are possible. In fact, because of the circular shape of one family of characteristics, 
the variables of the problem can be separated and solutions are regular. 

Regular inertial modes in a spherical shell therefore exist, but are these toroidal modes 
the only regular modes? we have no mathematical proof of it but numerical computations 
of the whole spectrum of eigenvalues including viscosity strongly suggest that this is 
indeed the case. The argument is as follows: regular modes in a spherical shell meeting 
stress-free boundary conditions have damping rates proportional to the viscosity which 
will turn out to be very small compared to those of singular modes which, as we shall 
see, develop shear layers. In a plot of the eigenvalues in the complex plane, regular modes 
will pop out when the viscosity is sufficiently low as can be seen in the context of gravity 



modes in Rieutord & Noui (199£ ). Computations for different m's show that only one 



eigenvalue popped out and that is the one of the toroidal mode. 

2.6. A summary of the results on the inviscid problem 

Before jumping into the question of how inertial modes of a spherical shell behave when a 
slight amount of viscosity is included, it is certainly useful to summarize the main results 
obtained in the foregoing sections on the inviscid problem. 



We have seen in §2.2 that the nature (regular or singular) of eigenmodes is, with the 
exception of toroid al m odes, determined by the dynamics of characteristics. The study 
of this dynamics (§ (2.2| ) revealed the generic property that characteristics converge to 
attractors made of a periodic orbit which exist in some frequency band. These attractors 
are also characterized by their length (i.e. the number of reflexions) which influence the 
rate at which characteristics converge to them; this rate is given by a negative Lyapunov 
exponent. When the frequency of the attractor is close to sin(p7r/2g) where p and q belong 
to a finite set of integers determined by the size of the inner core, the corresponding 
attractors are very long and weakly attractive. These frequencies are the ones for which 
the shadow of the inner core follows a periodic orbit (see figure ||); they will prove to be 
important in the determination of the asymptotic spectrum of inertial modes when the 
viscosity vanishes. 

The focusing of energy by attractors is not the only source of singularity: we have 
shown that near the critical latitude of the inner boundary a milder singularity will 
develop in general. This singularity will prove to be relevant in the viscous case when 
shear layers associated with attractors are inhibited. 

Finally, we found that some regular solutions still exist. They are purely toroidal modes 
and we surmize that they are the only true eigenmodes of a rotating fluid in a spherical 
shell. From the mathematical point of view, the point spectrum of the Poincare operator 
(i.e. eigenvalues associated with square-integrable functions) is almost empty. 



3. The solutions with viscosity 

3.1. General results 

When viscosity is included the equations are elliptic and the problem is well-posed; hence, 
the solutions are smooth C°°-functions which can be computed numerically. We shall not 
describe the method used and refer the reader to Rieutord & Valdcttaro (1997). We just 
recall that we solve the eigenvalue problem (A is the complex eigenvalue) 
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Xu + e z x u = -VP + E\7 2 u 
V • u = 



(3.1) 



with stress-free boundary conditions to eliminate Ekman boundary layers; E — v /2VLR 2 
is the Ekman number, v being the kinematic viscosity. 

The main result, coming from numerical solutions of this problem, is that the amplitude 
of the modes is concentrated along paths of characteristics drawn by attractors. However, 



as found by Dintrans et al. (1999) and Rieutord & Valdettaro (1997), the attractor 



appears in the viscous solutions only when the Ekman number is low enough. This 
critical Ekman number, below which the mode seems to reach an asymptotic shape, 
depends on the length of the attractor; short (and simple) attractors appear at higher 
viscosities than long (and complex) attractors which may never appear within the range 
of physically relevant Ekman numbers (E>10~ 1& ). 

In order to investigate the properties of viscous solutions associated with attractors, 
we shall focus on a few simple ones which appear at reasonable Ekman numbers (i.e. 
larger than 1CV 9 ). Some are the ones displayed in figure [?] plus two others located in 
the 0.6-0.625 frequency band, one of which was considered by [Israeli (1972 ) using a thin 
shell. 

A plot of the eigenmodes associated with these four attractors is shown in figure fL2[ 
This figure displays the kinetic energy of the modes in a meridional section of the shell. 
As expected the kinetic energy focuses around the attractors which we overplot on each 
diagram; however, this is not systematic as shown by figure p~2) c: there, the kinetic energy 
concentrates along a characteristic path starting at the critical latitude rather than along 
the (only) existing attractor. We understand this situation as the consequence of the 
location of the attractor: one of its segments is indeed almost tangential to the outer 
sphere which therefore inhibits the development of the shear layer. By computing the 
distance between the boundary and the attractor, we estimate that a £ ,1//4 -shear layer is 
inhibited by the boundary if the Ekman number is larger than 10 -11 . Such low Ekman 
numbers are out of reach numerically at the moment. This example illustrates the point 



mentioned in § 2.2.4 : very long attractors will appear at extremely low Ekman numbers. 
If we consider that lowest Ekman numbers are those of stars (~ 10~ 18 ) or the Earth's 
core (~ 10 -16 ) and if we use the same scaling as above for shear layers, then we can 
conclude (from figure ^) that attractors longer than ~ 100 will never appear in physical 
systems. 

We therefore see that, although the singularity associated with an attractor is stronger 
than that of the critical latitude, this latter singularity may show up if, for some reason, 
the build up of shear layers around the attractor is inhibited. We surmise that 'long' 
attractors, will dominate relative to the critical latitude singularity only at very low 
Ekman numbers. 'Short' attractors may therefore appear more easily as in figure [l2| a 
while still showing, weakly, the critical latitude singularity. 

Another surprising feature of the rays (i.e. shear layers) lying along an attractor is 
that the maximum energy density is not always centered on the attractor (figure |l2| b or 
|l2|d). We discuss this point below. 

To make some progress in the understanding of this complex behaviour we shall inves- 
tigate in more detail the structure of shear layers lying near the attractors. 



3.2. Structure of shear layers 
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Mode 0* j]=0.35 L=400 Nr=18D E=2.0xlo" 



Mode 0* jj=0.35 L=6B0 Nr=33D £ = 1.0x10" 



Figure 12. Distribution of kinetic energy in the meridional section of the shell for four different 
ax isvmmetric modes. These solu tions have been computed numerically using the same code as 
in Rieutord & Valdettaro (1997). On each panel, r is the damping rate, stress- free boundary 



conditions are used, L is the number of spherical harmonics and Nr is the number of grid points 
in the radial direction, (a) shows the mode associated with the equatorial attractor (in blue) of 
figure ^a while the polar attractor (in green) is only weakly visible, (b) Using a more damped 
mode with a slightly different frequency, we obt ain a case w here the polar attractor is fueled 
with energy. In (c) the attractor considered by Israeli (1972) should be fueled, but is in fact 
hampered by the boundary and the critical latitude singularity dominates the flow, (d) A very 
neat mode which concentrates along its attractor; in blue the corresponding equatorial attractor. 



3.2.1. Some numerical results 

As a preliminary step, we first compute the variations of the components of the velocity 
field along a line crossing a ray perpendicularly. Results are displayed in figure [l^. 
These profiles show that these internal shear layers have a rather complex structure 
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Figure 13. Left: the real and imaginary parts of V v in the cross-section of the attractor displayed 
in figure [12a. The segments with positive slope have been shown in cross section. The + sign 
overplottea on the curves represents the variations of iVn; the perfect matching with the curves 
of V v sho ws t hat t he phase quadrature between these components is well verified as expected 
from (3^) or ( 3.13 ). The dashed line shows the amplitude of the \v\ profile and the two crosses 
on the j/-axis indicate the position of the attractor. Right: same as on left but for a mode with 
more complex rays: it is a cut through the rays with negative slope of the mode of figure [l2|d; 
the cut starts near the critical latitude and is perpendicular to the rays. 



which looks like a plane inertial wave trapped in a 'potential well'. Each mode seems 
to be characterized by the number of nodes in the cross-section of its rays, just like a 
solution of a Sturm-Liouville problem. The analogy cannot, however, be pushed too far 
since the actual oscillations do not disappear outside the rays but continue with a very 
low amplitude (the well is leaking!). This is a consequence of the fact that the 'well' is 
not a local well but the result of a mapping made by the convergence of characteristics 
towards the attractor. We also surmise that since the convergence rate is not the same 
on each side of the attractorff], the 'potential well' is certainly not symmetric with respect 
to the attractor; we thus explain our finding that the maxima of kinetic energy density 
are not centered right on the attractor as shown by figures |l2|b,d or [l^. 

In the above view, the shear layers result from a balance of the focusing action of the 
mapping and the 'defocusing' action of viscosity; because the former action is global and 
the latter is local, the boundary layer analysis is difficult, if not impossible. The following 
analysis gives some general properties of these shear layers, properties which are actually 
observed numerically, but is not able to reproduce their detailed structure. 

3.2.2. Boundary layer analysis 

In order to describe the shear layers featuring the inertial modes of a spherical shell, 
it is convenient to project the equations in the characteristics' directions. 

f We mean here at some finite distance from the attractor; right on the attractor the conver- 
gence rate is given by the Lyapunov exponent. 
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From (2.7), we derive the expressions of unit vectors parallel (|j) or perpendicular (_L) 
to characteristics of positive (+) or negative (— ) slope: 



ej. = use s ± ae z , 



■ we z =F ae s 



Restricting ourselves to the case of a characteristic with positive slope, the components 
of the velocity in a meridional plane are: 

Vj| = ljVs + aV z , V± = -aV s + loV z 
We may now transform the equations of motions, written in cylindrical coordinates, 

dP 



XV S -V V = -—+EA'V 6 

XVp + V S = EA'V V 
dP 



XV, 



dz 



(3.2) 



into 



AVu - ujV v 



dP 

dx 



E V 



acuE . 



-V± 



XV V + uV\\ - aV± - EV 2 V V 
dP 



dy 



E V 



V± 



aioE . 



-V\\ 



(3.3) 



where A' = V 2 — 1/s 2 , and x and y are local coordinates respectively parallel and 
perpendicular to the characteristic. For the sake of simplicity, we may consider one such 
characteristic so that 



Thus 



x = az + cos, y = ujz — as, s = ljx — ay 



d d d d d d 

os ox oy oz ox oy 



(3.4) 



Mass conservation requires that 



dV \\ 
dx 



8V± 
dy 



u>V\\ — aV 



J 

s(x,y) 



± 







(3.5) 



As it was shown in Ricutord & Valdcttaro (1997), the inviscid balance along rays 
shows a dependence of the velocity and pressure fields with 1/y/s; we shall remove such 
a dependence from our equations by setting V = u/\fs and P — p/y/s. Hence, (3.3) and 
©yield 
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Figure 14. Shape of the Moore-Saffman function with m = — 1/3. 



dp ujp , 2 



OLLjE 



\u± + au v - 
du\\ du± uu\\ — au± 



oiloE 



dx dy 



= 



where now V 2 = + -jr-* 

as 2 4s 2 



2s 

d 2 _ a 2 

dz 2 dx 2 dy 2 4(ujx — oty) 2 



3 2 



(3.6) 



3.2.3. The inner E 1 ^ 3 -layer 

Searching for a boundary layer solution scaling with E x l 3 , we make the expansion 



(3.7) 



u|| = 4 + E 1 / 3 u\ + ■■■ 
u v = u$ + E^uf + ■ ■ 

u± = E 1 ' 3 ^ H 

p = E 1 / 3 p l H 

We shall use the scaled variable Y = yjE 1 ! 3 . We also recall t hat A = iuj + t, [lu, t) G H 2 
and that |r| <C From the second and third equations of (|3.6|) at zeroth order we get 



* ~ — 



and 



9pi 



"o - '"o ""o qy 

while from the combination of first order terms we get 



(3.8) 
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Figure 15. Motion of the virtual source after a reflection: the source A moves to A' after the 

reflection of the ray. 



-lot 



(3.9) 



dY 3 dx 

Making a last change of variables q — x/a and dropping the zero-index, we finally obtain 



,du u 



Q Y 3 



(3.10) 



which was first derived by Moore & Saffman (1969) for steady (vertical) shear layers. 

Moore and Saffman have shown that the solutions of ( 3.10| ) which can describe a 
detached shear-layer are self-similar solutions of the form: 



= q m H m [Y/q 



tV3 



(3.11) 



where the function H m is defined by 

H m (t) -- 



e~ lpt e 



p p-'^dp 



Since (3.11) describes detached shear layers, the function H m needs to vanish when 
t — > ±oo which is possible only if m < ( Moore fc Saffman 1969 ). The shape of this 
function is given in figure |TJ. 

To complete the description of the -| layer we need to determin e the index m of the 



Moore and Saffman function. For this purpose, we first note from ( |3.11 ) that the width 
of the layer is singular at the origin of the ir-axis. This origin can therefore be considered 
as a virtual source of the ray which obviously lies outside the fluid's container. Let us 
therefore consider a segment of a mode around an attractor. Let us orient the x-axis in 
the direction of contraction of the map and call xo the abscisse of the first point of this 
segment (point B in figure 111) . At B the width of the layer is proportional to x^ 3 while 
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the amplitude of u v is proportional to x™H m (0). At C and before reflection, the width 
is (xo +t} 1 ^ and the amplitude is (xo + £) m H m (0); after reflection, the mapping changes 
the scale by a factor K; therefore the width is now (xo + £) X I^K and the amplitude 
(xq + £) m H m (0) / K . The width is just as though the virtual source A' were at a distance 
(xq + £)K 3 of C while the amplitude would imply a distance (xq + £)K~ x l m \ since the 



segment starting at C must also be a solution of the form (3.11), the new virtual source 
must be at the same position for both the amplitude and width; therefore, we need to 
have 



m = ~ (3.12) 



This index was also found by Stewartson (1972a) on the argument that it is the only 
one for which the flux 

V ll dY 

■oo 

is conserved along the ray, which means that the ^ layer does no pumping. 

From the property that H m (t) ~ t 3m as t — > oo, we see that the solution (3.11) is 
independent of x far from the layer and decreases as 1/y. 

3.2.4. The outer E l ■> /A -layer 

Since some mode s show a clear scaling of their ray with the \ exponent (see Rieutord 



fe Valdettaro 1997), we also briefly discuss this case. 

As the j -layer is also much larger than the Ekman layer, the ip- and || -components are 
in quadrature, i.e. , 

V v = iV\\ (3.13) 
These components also verify, from mass conservation and inviscid balance, 
V v =F(y/EV*)/^. 

We cannot say much more about the ^-layer, except that the expansions to higher 
orders including viscous terms lead to a differential equation for F which is not closed, 
some additional functions missing their differential equation. If these extra and undeter- 
mined functions are set to zero, F obeys, as expected, a fourth order differential equation 
whose solutions do not agree (in general) with the numerical results. We think that this 
is due to the action of the mapping which is obviously missing in the local analysis; 
unfortunately, we have not yet found a way to include it. 

3.2.5. Wave packet kinematics 

In order to give a more physical understanding of the behaviour of shear layers as 
viscosity is reduced, we propose to consider a wave packet traveling around an attractor. 

Let us suppose that the Ekman number is very small but finite. When traveling along 
an attractor a wave packet is damped by viscosity but its reflections on the boundaries 
enhance it if the direction of propagation is such that the map is contracting (A < 0); 
this equilibrium may be written as: 

e -£„*;Un =e .vA (3.i4) 

where N is the length of the attractor, A its Lyapunov exponent^, t n the time elapsed 

f recall that according to the definition of the Lyanunov exponent, e A is the mean dilation 
rate (of an interval 8(j>) per bounce. 
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on the n'^-segment and k n the wavenumber of the n t/l -segment. Noting that on each 
segment the group velocity is almost constant and reads: 



we may transform (3.14) into 



A = - 



nVT 



(3.15) 



where we have introduced the length £ n of each segment of the attractor; we also used 
the fact that k s — k\/\ — ui 2 so that t n = £ n k n /\/l — ui 2 and all quantities are now 
dimensionless. Since k n = C n fc n -i, C n > 1 being the contraction coefficient of the n th 



reflection, we may rewrite (3.15) 



A 



Ekf 



N 



where 



N vT 



F{u) 



Y,^Y[Cf = -EkfF(u 



(3.16) 



n=l i=2 



1 N n 



NVl - LO 2 



n=l i=2 



is a purely geometrical quantity describing the path of characteristics associated with the 
attractor at the frequency uj. The expression (3.16) expresses through ( |3.14 ) the strict 
periodicity of the amplitude of the velocity field along an attractor when viscosity is 
small but finite. 

A similar relation may be derived if we now express that the scale of a wave packet must 
be the same after one cycle along the attractor. In a purely diffusive (viscous) process, 
the scale of a structure, initially being 'a' grows like Va 2 + vt with time; therefore the 
relation between the layer's width after one propagation and one reflection is 



a n = D n \ a 2 n _ x + vt r 



(3.17) 



Here D n is the dilation coefficient of the n th -reflection (D n = 1/C„ < 1). Along one cycle 
with N reflections, we have 



i N 



Using the fact that a^v+i = ai, (3.17) leads to 



A = -J-Vlnfl + % > ) 



imposing 0<(7<l/3, we finally obtain 



A = - 



E 



2NVl 



P k 



Noting that a n ~ A n = 27r/fc ra , we recover (3.15) except for a constant factor. 

The two derivations of the Lyapunov exponent through this schematic model show that 



(3.18) 



Now, if we let the width of rays a n scale with E a , we find that vt n la? n ~ E 



1-3(7. 
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Figure 16. Left: The inverse of the length N of all the attractors with a length less than 100 
for a spherical shell with n = 0.35 are represented with a (*) denoting the frequency (or critical 
latitude) where A = and with a line segment showing the interval of existence. By showing 
the projection of the * on the x-axis, we try to give an idea of what would be the asymptotic 
spectrum neglecting very long attractors. The two vertical lines show the position of tv/8 and 
7r/6 which correspond to the periodic orbits of the shadow for n = 0.35. Right: A blow-up of 
the region around 7r/6; note the lengthening of the attractors as this value is approached. 



the width of shear layers lying along a periodic attractor and scaling with E a , should be 
such that 

1 

CT< 3 

We therefore see that the ^-exponent is a limit case. In fact this inequality shows that 
'naked' ^-layers cannot exist and should be embeded in thicker layers; this seems to be 
the case indeed, at least for all the modes which we investigated in detail: they usually 
show (7 — 1/4. 

3.3. The asymptotic spectrum 
The foregoing calculations show one important result: As viscosity tends to zero and 



since k\ cx E~~ a , from ( 3.16 ) we may conclude that the Lyapunov exponent of attractors 
must vanish as viscosity vanishes following the law A cx E 1 ~ 3a . It therefore turns out 
that cigcnfrcqucncics will converge towards the roots uii of the equation A(u>) = which 
therefore describe the asymptotic spectrum of inertial modes in a spherical shell. From 
the fact that only a finite number of attractors exist at a given frequency, we deduce that 
the spectrum cannot be dense in [0,1] contrary to the case of the full sphere. However, 
this spectrum has some accumulation points ui a which are due to the existence of neutral 



(A = 0) periodic orbits with frequencies sinp7r/2g (see §2.2.2). Indeed, in the neighbour- 
hood of such points we may find attractors which are longer and longer, the closer they 
are to u) a . However, the number of accumulation points is finite and given by the number 
of pairs (p, q) possible for periodic orbits of the shadow and when r\ > 1/ V2 only three 
such points (u> = 0, 1/ \/2, 1) exist. In figure [l6| we clearly see the accumulation points cor- 
responding to sin(7r/4) and sin(7r/6)(j]. When 77 — ► 0, the number of accumulation points 
gets larger and larger, as more and more rationals are added to the set of accumulation 



f The case sin(7r/8) is not as clear for it needs a much higher frequency resolution since the 
shadow almost fills the whole volume as it/8 ~ arcsin(0.35). 
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Figure 17. Distribution of the 
least-damped eigenvalues in the com- 
plex plane when E — 10 -8 with resolution 
Lmax=700 and Nr=270. The dotted line 
shows 7r/6 while hatched bands indicate 
the positions of three simple attractors; 
n = 0.35. 



Figure 18. Asymptotic behaviour of the 
eigenvalue associated with the eigenmode 
plotted in figure |l2^ ,. The dashed line repre- 
sent so; — Uas where u>aa. is in fact u>\ given 
by ( |B 2| ) (cf appendix [b]), while the dotted 
line is for the damping rate. The solid line 
represents the 'theoretical' law E 1 ^ 2 . 



points, a situation in accordance with the fact that at rj — the spectrum is dense in 
[0,1]- 

Let us also underline the fact that when E = 0, eigenvalues disappear since solutions 
of the equations are no longer square-integrable; this is also true for frequencies of at- 
tractors such that A = 0, since attractors still focus the energy (but algebraically, not 
exponentially) . 

For a given upper bound of the damping rate, eigenvalues will be packed around 
the roots w% and around the allowed frequencies of the set sin(p7r/2g). In figure [17], we 
computed the distribution of least-damped eigenvalues when E = 10~ 8 , i.e. for 140 
evenly spaced frequencies between and 1/V2, we computed the eight least-damped 
modes. This figure offers a glimpse at the asymptotic distribution of eigenvalues in the 
complex plane: we clearly see three main bands[f] of attractors where modes are more 
damped and the two frequencies sin(7r/6) and sin(7r/4) where least-damped modes tend 



to accumulate; the sin(7r/4) case is conspicuous. Note also the similarity with figure 16 
where the three bands made by the aforementioned attractors are also clearly visible. 

From the asymptotic behaviour ( 2.16| ) of the Lyapunov exponent in the vicinity of the 
roots uji, we can derive that 



uj = lu, + aE 2 - 6a + ■■■ 

while an order of magnitude evaluation of the ratio of dissipation to kinetic energy yields 
the asymptotic law of the damping rate r 

r = -bE 1 - 2a + ■■■ 

f They are [0.3959, 0.4162], [0.5290, 0.5554] and [0.6 ,0.6266]: the first and third are illustrated 



by attractors in figure h2t the second is illustrated in Rieutord et al. (200C); for all n = 0.35 
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This asymptotic behaviour of eigenvalues is best illustrated in figure |l8| where 
a = 1/4. Such a mode is the least-damped one in its frequency range and is therefore 
not perturbed by other eigenvalues. This is likely the reason why it shows its asymptotic 
regime at rather 'high' Ekman numbers. 

Concerning the eigenmodes, it is worth noting that (2.16) and (B9) imply that the 
spacing of adjacent rays scales like E X I A (if a = 1/4); therefore, the distance between two 
rays of an attractor remains the same when it is rescaled by E 1 ^ and one may conclude 
that each mode in the form of a viscous attractor keeps a self-similar structure as the 
Ekman number vanishes. 



4. Discussion 

Ending this paper, we think that the asymptotic behaviour of inertial modes in a 
spherical shell when the Ekman number vanishes can be anticipated even if some points 
remain in the shadows. 

We have seen at the beginning of the paper that the trajectories of characteristics in 
general converge towards an attractor; exceptions are when the sphere is full or the inner 
core is small enough to let a finite number of periodic orbits remain (which are associated 
with critical latitudes commensurable with tt) . Leaving the full sphere for which analytical 



solutions exist since Bryan (1889), the generic behaviour of characteristics is that they 
converge towards an attractor which is a periodic orbit residing in some frequency band. 

The knowledge of the characteristic trajectories can be used immediately in two- 
dimensional problems in order to construct a solution of the inviscid problem. This 
solution contains an arbitrary function which needs to be specified on some fundamental 
interval(s). This makes the eigenvalues always infinitely degenerate. In three dimensions, 
the trajectories cannot be used so efficiently but their convergence towards an attractor 
can be used to show the divergence of the velocity field at zero viscosity. In two dimen- 
sions, this divergence allowed us to prove the non-square-integrability of velocity fields 
associated with attractors, implying the absence of eigenvalues in a large fraction of the 
frequency band [0,2O]. 

Beside the singularities generated by attractors, we also shed new light on the singu- 
larity arising at the critical l atitude of the inner shell. We thus generalized the result of 



Btcwartson fc Rickard (1969| ) that the velocity field diverges as the inverse of the square 
root of the distance to the characteristic grazing the inner shell; this singularity also 
makes the velocity field not square-integrable. 

Among all these singular solutions, a small set of regular modes 'survive': they are 
purely toroidal modes which, thanks to a velocity field which has no radial component, 
do not suffer the constraints imposed by characteristics paths. Numerical computations 
of the whole spectrum give a strong evidence that these modes are the only regular ones. 

When viscosity is included, all the aforementioned singularities appear in the form of 
shear layers. In the asymptotic regime, we therefore expect that attractors will feature 
the viscous solutions. However, this asymptotic regime may be reached at extremely 
low values of the Ekman number, some of which may not even be relevant astrophysi- 
cally or geophysically. We may therefore face some intermediate regime where the milder 
singularity at critical latitude plays an important part in featuring the amplitude of a 
mode. 

Our numerical investigations of the structure of shear layers which are generated by the 
different singularities revealed a rather complex structure of nested layers scaling with E a , 
< g < 1/3 where the value a =1/4 seems to be favoured. In some simple cases, where 
the velocity field shows no node in the transverse direction of a shear layer, numerical 



Inertial waves in a rotating spherical shell: attractors and asymptotic spectrum 33 

results indeed show a scaling with E 1 ^. Our boundary layer analysis demonstrated that 
all these internal layers should contain an inner a = ^-layer which is similar to the 
vertical Stewartson layers; however, the mapping made by characteristics influences the 
scales larger than E 1 / 3 and therefore makes a local analysis insufficient to determine the 
structure of the outer parts of the layers. 

The upper bound a < 1/3 has been derived using a heuristic model of an inertial wave 
packet traveling along an attractor; from this model, we also showed that the asymptotic 
spectrum of eigenvalues can be derived, once the structure of the shear layer can be 
computed. An important result of this analysis is that the limits of eigenfrequcncies as 
E — > 0, do not form a dense set in [0, 20] contrary to the case of the full sphere. 

The asymptotic behaviour of inertial modes and their associated eigenvalues is there- 
fore slowly becoming clearer: As the Ekman number decreases, more and more eigen- 
modes are concentrated along the attractors associated with their frequency; once this 
asymptotic regime is reached the frequency of the mode changes slowly with viscosity 
so that the Lyapunov exponent of the attractor decreases (in absolute value) and con- 
verges toward the frequency where this exponent is zero. We are thus in the position of 
describing regimes with extremely low values of the Ekman numbers that are relevant in 
astrophysics (E = 10~ 18 for a radiative zone of a star) or geophysics (E = 10~ 15 for the 
liquid core of the Earth) and which are way out of reach numerically. 

However, some points remain in the shadows, the most challenging one being the 
structure of the shear layers. As we observed, this structure builds up from a large scale 
phenomenon which is represented by the mapping of characteristics and a small-scale 
one which is diffusion. To the best of our knowledge, such a problem has never been 
investigated in the past. Since the three-dimensional case is much more involved than 
the two dimensional case, because of the intrusion of Riemann functions, we think that 
this latter case should be investigated first; this will be the subject of future work. 

The foregoing results were derived from the analysis of inertial waves of a fluid con- 
tained in a spherical shell but it is clear that they are of wider generality. They can 
be easily generalized to any container of the same topology, like ellipsoidal shells, or be 
used qualitatively for any kind of container. In the case of ellipsoidal shells, even axial 
symmetry can be relaxed: since constraints imposed by characteristic surfaces are in a 
meridional plane, no small-scale should appear in the ^-direction. Attractors are robust 
structures and only the longest ones are sensitive to small modifications of the shape of 
the boundaries; however, as they would transform into other long attractors, the final 
solution would not be much affected. This kind of 'structural stability' is important when 
applying these results to real objects like the core of the Earth which is obviously not a 



perfect spherical shell (see the discussion in Rieutord 2000a) 



Our results also naturally extend to all systems governed by a spatially hyperbolic 



equations. Hence, one will find similar properties for gravity modes (Maas & Lam 1995 



Rieutord fc Noui 1999 ) or hydromagnetic modes ( Malkus 1967 ). 

Now the next question raised by our results concern their implications when the modes 
are of finite amplitude. These may concern the development of the elliptic instability since 



this instability is precisely an instability of inertial modes (see Rieutord 2000a ) or may 



affect the transport properties of the fluid which arc much enhanced around attractors 



( [Maas et al. 1997J pintrans et al. 1999j ; [Dauxois fc Young 1999|) 



In the same context, one may wonder whether the attractors can be studied exper- 
imentally. At the moment, only one attractor has been detected experimentally, using 



gravity waves of a stably stratified fluid (Maas et al. 1997) or inertial waves (Maas 2000) 



The main obstacle for detecting attractors with experiments is the rather large value of 
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the Ekman number of experiments. Numerical calculations have indeed shown that this 
number should not exceed a few 10 -8 . Using water in a spherical shell with a radius of 
20 cm demands a rotational speed of 12000 rpm which is difficult to achieve and raises 
experimental problems. Also attractors should not be confused with other phenomena 
which emphasize characteristics paths like the critical latitude singularity or a forced 
perturbation like a discontinuity in velocity forced by boundary conditions (e.g. the split 
disk case). 

Finally, these new features of inertial modes may also have some interesting conse- 
quences in astrophysics. It is now well-known that rapidly rotating neutron stars can 
lose a substantial amount of angular momentum when some inertial modes become un- 



stable beca use of a coupling with gravitational radiation (see |Andersson 1998] ; |Lindblom 



et al. 1998 ). This instability therefore controls the rotation speed limit of neutron stars 
and this limit would be the higher, the more damped are inertial modes. Stars with a 
core or density jump will therefore be more stable than others, a fac t which may be u sed 



to give new constraints on the state of matter inside neutron stars ( Rieutord 2000b ) 
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Appendix A. arcsin(\/3/7) 

In this appendix, we show in which cases sin(pir/q) is the square root of a rational. 
This establishes that y/Z/7 is not the sinus of any rational fraction of n. The proof may 
exist in the mathematical literature, but we have not been able to locate it; we therefore 
propose here a simple proof of the result. 

Let p/q be a rational number, p and q being coprimes, and let us suppose that 
sin.(pn/q) = \J a/b, where a/b is a rational number. Then cos2pn/q = 1 — 2sin 2 pn/q 
will be a rational number. 

Let us take first the case where q is prime. Then exp(2«p7r / q) is a q th root of unity, 
and as such is a solution of: 

X 11 - 1 + X q - 2 + ■ ■ ■ + 1 = 0, (Al) 
provided q > 1. But one has also: 

9-1 

Xi- 1 + x q - 2 + ■■■ + != Yl (X- e l2m7r /«) , (A 2) 

m—l 

since the exp(i2m7r/g), m = 1, . . . , q — 1 are roots of the polynomial (called cyclotomic 
polynomial). 

Since exp(i2pir / q) is root of this polynomial, so is exp(— i2pir/q), and they are different 
if q > 2. The product ( [A 2| ) thus contains 

(X - cxp(i2 P Tr/q)){X - cxp(-i2 P TT/q)) = X 2 - 2X cos(2pw/q) + 1 (A3) 

Therefore, since cos(2p7r/g) is rational, X 2 — 2X cos(2pTr/q) + 1 is a rational polynomial, 
which divides the cyclotomic polynomial X q ~ x + X q ~ 2 + . . . + 1. But Gauss proved (see 
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Jacobson (1985 ) p. 272) that this polynomial is irreducible in the field of rationals if q is 
prime. Therefore the degree of this polynomial, which is q — 1, cannot exceed two, since 
otherwise it would be reducible. So q = 1, 2, 3 are the only possibilities. 

If q is not prime, then the polynomial X^ 1 + X q ~ 2 + . . . + 1 is not irreducible. But 
exp(2ipir / q) is still a q th root of unity. If exp(2ipmir / q) ^ 1 for all m < q, one says 
that exp(2ip7r / q) is a primitive q th root of unity. exp(i2p7r/q) is a primitive q th root of 
unity since p is prime to q (Hardy & Wright (1975)). Then obviously exp(— i1pi:jq) is 
a primitive q th root, different from exp(i2pTT / q) since q > 2. Therefore the polynomial 
whose roots are all the primitive q th roots of unity contains the factor (A3). But this 
polynomial is ir reducible in Q (cyclotomic polynomial) (see Hardy & Wright (1975), 
Jacobson jlggf )). So if the degree of the polynomial is greater than two, there is a 
contradiction. The degree of the polynomial is equal to the Euler function faq) which is 
the number of positive integers not greater than and prime to q. One has faq) > 2 for 
q > 6. 

Therefore the only possible values for q are q = 1,2,3,4,6. By direct verification, one 
sees that \j2>j7 is not sin(p7r jq) for one of these q. 



Appendix B. Lyapunov exponent of two attractors 

B.l. Some example of Lyapunov exponents 

B.l.l. Parameters of the equatorial attractor 

The equatorial attractor drawn in figure 0a exists for all A's such that the following 
equation has a root for fa € [0, A[, the latitude of the point of reflection on the inner 
sphere, 

6A = arccos(?ycos(A + fa)) + arccos(?7 cos(A — fa)) 
as can be derived from the relations 



fa + fa = 2A 

fa + fa = -2A 

cos(fa + A) = r]cos(fa + A) 

cos(fa — A) = r)cos(fa — A) 




The frequency band [a>i,c<j2[ where this attractor exists is such that: 



V 



and 



U>2 



sinA2, with cos(6A2 — arccosry) = 77cos(2Aa) (B2) 



If i] = 0.35, we find w x = 0.403112887 and u 2 = 0.412474677. 
For this orbit the Lyapunov exponent is given by 



A(w) 



1 



InCV) 



(B3) 



where C(u>) is the contraction coefficient (this orbit is attracting when it is followed in 
the trigonometric sense) . Using the colatitudes of the reflection points, it can be obtained 
with 



C — CxC^C^C^ 



A) sin(0 2 + A) sin(0 3 — A) sin(0 4 + A) 



+ A) sm(fa — A) sin(0 3 + A) sin(fa — A) 



If we use the fact that fa + fa = 2A and fa + fa = 



—2 A, we finally get 
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AM 



1 



In 



6 3 - A) sin(</> 4 + A) 



sin 



(B4) 



63 + 5A) sin(04 — A) 

Such a relation can also be obtained by differentiation of the formulae relating the 
angles 0i, 03 and 04 in the third and fourth relations of (Bl); this yields: 



sin 



t> 3 - A) sin(0 4 + A) 



t>3 sin(0 3 + 5A) sin(0 4 - A) 



B.1.2. The associated polar attractor 

Reflection points of this attractor (left in figure fjja) are related by: 



4>\ + 4>a = 2A 

02 + 03 = 2A 

03 + 04 = -2A 
COS(02 + A) = 77COs(05 + A) 

— cos(0i + A) = 7/cos(05 — A) 



which implies solving 



(B5) 



(Bl 



8A = arccos(?ycos(05 + A)) + arccos(— 77 cos(A — 05)) (B7) 
for A < 05 < 7r/2. The bounds of the interval of existence are [Ai, A2] such that 

cos(8Ai — arccos(— 77)) = ?7cos2Ai 

cos 4A2 + 77 sin A2 = 

Here A(Ai) = -00 , A(A 2 ) = 0. When 77 = 0.35 we find uji = 0.395915 and u 2 = 0.416185. 

In a similar way as we derived (B4), we find that for the polar orbit, the Lyapunov 
exponent reads 



A( W ) 



hi 



sin(A + 0i ) sin(A + 5 



sin(7A — 0i) sin(A — 05) 



(Bl 



B.2. A in the neighbourhood of the point where A = 
B.2.1. General result 

In this subsection we prove that generically near a critical latitude Ao which displays 
a periodic orbit having A = 0, the Lyapunov exponent behaves like the square root of 
the variation of the critical latitude. Actually, this is a general result valid for a one- 
dimensional mapping that depends on one control parameter (A in our case) near a 
tangent bifurcation point, provided the mapping is sufficiently smooth there. 

Let /(0, A) be the mapping (for cj), in our case it is f N (<fi) with / defined in ( 2.11 ) and 



is the period of the orbit) depending on the parameter A. Let Ao be the value of A for 
which A = 0, and (j)Q the fixed point of the corresponding map (/(0o, Ao) = 0o). Ao is a 
bifurcation point, as illustrated in figure ^9[ and therefore <9//<90|0 o .a o = 1. 

Let 0i be a fixed point for the mapping at Ai near Aq: 0i = /(0i, Ai). Let us define 



fa 



d l+j f 



SX — Ai — Ao 



00,-^0 
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ro q>T $ 

Figure 19. Generic behaviour of the mapping near a tangent bifurcation point (</>o, Ao). 

A Taylor expansion of / around (fa, Ao) up to order 2 yields: 

Ai) = f(fa, Ao) + /io<ty + foiSX + \h^ 2 + l.fo2SX 2 + fn6\6<f> + ... 
Since f(fa, Ai) = fa, f(fa, A ) = fa, and f w = 1, we get: 

./• 20 <5</> 2 + + 2/oi<5A + f a2 6X 2 = 

Solving this equation for S(f> gives: 

§(f> = -/nSA ± ^(/^ - fwf^SX 2 - 2f 01 f 2Q S\ 

/20 

for the coordinates of the two fixed points at Ai, <pi and 0f . 
To leading order in SX we obtain: 



5<f> ~ ±J -2^VS\ 

V /20 



(B9) 



That is, the displacement of the fixed point 5(f) behaves like the square root of the 
variation of the control parameter A. 

We compute finally the Lyapunov exponent associated with the mapping at (fa, Ai): 



A = In 


d.f 




~ In I/10 + /20<5</> 


+ fuSX + ... 


~ In 




d<f> 


<t>i, Ai 









i± f Jz^iVsx 
V /20 



20 \ 



-2/< 



01 



/20 



The positive value corresponds to the repulsive fixed point </>f = fa + y — 2j^VdX, 
while the negative value is associated with the attractive fixed point 

fa 4 - = fa-^/-2^V6X. 

Therefore it follows that if A(u>q) = for an attractor, then, 



' /20 



A(ui) ~ —Ay/\uj — LO \ 



in the neighbourhood of u>o- 
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B.2.2. Example of the equatorial attractor 



Let us differentiate relation (B3) 

dk 
dm 



dC 



4C cos A dX 



where C is given in (B 5) for instance. We introduce now the new variables a = 4>a and 
= 03 so that C reads 



C 



A) sin(a + A) 



sin(5A + <f>) sin(a — A) 



We now evaluate ^ at the point where A = 0; let us call Aq this point, it follows that: 



dC 

— (A ) = 2a' cot A - 2(2 + 0') cot 3A 
dX 

where a' = da/dX and 0' = d<f)/dX. We used the fact that Ao verifies 

cos 3A = r\ cos Aq 



(BIO) 



Using now ( B 1 ) , we have 

cos(5A + 0) — -q cos(A + a) 
cos(A — </>)= 7] cos(A — a) 
which we differentiate with respect to A 



(5 + (j)') sin(5A + </>)= rj{\ + a') sin(A + a) 
(1 - 4>') sin(A -4>) = 77(1 - a') sin(A - a) 

and solve that system for a' and <f>' . Its determinant is 



(Bll) 



(B12) 



which yields 
in the vicinity of Ao- Setting 

it turns out that 



A = ?/sin(A - a) sin(5A + 0)(1 - e~ 4A ) 
A ~ 4A?7 sin A sin 3A 



— — , and a = — — 
A A 



iV (A o ) = 277sinA o /(Ao)^O 

A^ Q (A ) = 2sin3A /(A ) ^0 

with /(A) = 77 sin A — 3sin3A. This shows that a' and <j)' tend to infinity in Ao- After 
substitution it turns out that 



Hence 



/(A )= -(3 + v)Vl Z7 r } 

dC cosAo , . , . 4(77 + 3) 

dX (Ao) = — H{Tl) Wlth H{V) = -WT^ 



Therefore 
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f"(if 2 ) f-«(t> 3 ) t> 2 ^ *i * 



Figure 20. Generic illustration of the mapping /" in the neighbourhood of an attractive fixed 
point 0i. (j>2 is the point of discontinuity nearest to 0i. The interval / = [/ _n (02), / _n (03)[, 
after application of the mapping enters the basin of attraction of 0i. 



dA 

dui 



Hijj) 
4A 



and finally 

At u>i we thus have 



A 



(W - LUl) 1 / 2 



vV^~~v V (v + 2) 2 



1/2 



If 77 = 0.35, then G(tj) = 4.8184134. 



Appendix C. The number of attractors at a given frequency 

We show here that the number of periodic orbits for a given frequency is bounded by 
the number of discontinuous points. 

To demonstrate this point, let us suppose the mapping has p discontinuities. This 
means that f k will have kp discontinuities, since the image of any discontinuous point is 
a discontinuous point. This means that f k is C°° on kp interval. Each interval is bounded 
by two discontinuous points. 

When an attractor of period n appears, it means that the orbit bounces n times on the 
outer shell. Therefore, n intervals in the graph of /" will cross the straight line y = x, 
where they will be locked in subsequent iterations of /". 

We shall denote by 4>i the attractive point closest to the end of the interval, and by 
4>2 the nearest discontinuous point of /", which bounds the locked interval; thus, [4>i, §1 [ 
belongs to the basin of attraction of Now, <j>i is a fixed point of /" and therefore 
of any iteration of /" or f~ n . On the contrary, <t> 2 is not. It is a point of discontinuity 
of therefore f~ n ((f>2) is a point of discontinuity of f 2n . Thus, f~ n ([4>i, 02 [) belongs 
to the basin of attraction of f 2n ; since, in general, f~ n {(j)i) ^ 02, f~ n {[<fii, <j>2[) is not a 
continuous interval and we see that the basin of attraction of f 2n contains at least two 
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intervals: [</>i, </>2[ and some other interval in the neighbourhood of f~ n (4>2) (for instance 
I = ]f~ n {4>v)i /~ n (<fe)] in figure |20| ). Therefore at the stage 2n, one of the new intervals 
created falls into the basin. This is true near all the n attracting points of f n , so n 
additional intervals fall into the basin at the stage 2n. After n other iterations, f~ 2n ((j>2) 
will be a new point of discontinuity of / 3n , not present in /" and f 2n , and the same 
argument shows that n additional intervals at least fall into the basin at stage 3n. 

This indicates that for each n iteration, n additional intervals (at least) are in the 
basin of attraction of the attractor of period n. Let us therefore consider a case with 
two attractors of period n\ and 712; after n\ iterations we get n^p intervals bounded by 
discontinuities but the attractor has captured n\ intervals; therefore outside the basin of 
attraction of the first attractor, we have, for the ni-iterate, ri\(j> — 1) 'free' intervals at 
most; if the mapping is iterated ri\ni times, we have n\ni(j>— 1) free intervals. These free 
intervals contain the basin of attraction of the second attractor; but after ri\rii iterations 
the second attractor has captured n\n2 intervals, therefore only n\n2{p — 1) — n\n2 = 
n\n2{p — 2) are really free. Following this reasoning for a third attractor, we would get 
n\n2n^{p— 3) free intervals left. Thus not more than p attractors can exist simultaneously. 

It is interesting to note that numerically the number of attractors was always found 
smaller or equal to p/2. We note that the argument above is actually an upper bound. 
In practice, the preimage of [a^j^i] can include other points of discontinuity already 
existing, and more than n intervals can fall in the basin after each n iterations of the 
mapping. 

The above argument is valid provided there is only one attractive point by interval 
between two discontinuous points. Actually, we never found through extensive numerical 
simulations of the mapping, a case where several attracting orbits coexist on the same 
interval (apart for the values of the frequency where families of neutral orbits exist). 
Actually, such a case would be related to a period-doubling bifurcation which cannot exist 
in our system since one cannot cross the straight line y = x with a negative derivative. 
Even if several orbits can coexist on a given interval, this will happen at a finite iteration 
of the mapping and the total number of periodic orbits will remain finite, since the 
argument above shows that the number of such intervals is finite. 
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